Free Energy Functional for Nonequilibrium Systems : An Exactly Solvable Case 
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We consider the steady state of an open system in which there is a flux of matter between two 
reservoirs at different chemical potentials. For a large system of size N, the probability of any 
macroscopic density profile p(x) is exp[—NT'({p})]; T thus generalizes to nonequilibrium systems 
the notion of free energy density for equilibrium systems. Our exact expression for T is a nonlo- 
cal functional of p, which yields the macroscopically long range correlations in the nonequilibrium 
steady state previously predicted by fluctuating hydrodynamics and observed experimentally. 
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The extension of the central object of equilibrium sta- 
tistical mechanics, entropy or free energy, to nonequilib- 
rium systems in which there is a transport of matter or 
energy has been the holy grail of nonequilibrium statis- 
tical mechanics since the time of Boltzmann. An impor- 
tant step in that direction was taken by Onsager and 
Machlup (!]] for linear deviations from equilibrium, and 
there have been many further extensions to time depen- 
dent evolutions starting and remaining in local equilib- 
rium [0-0]. The extension to stationary nonequilibrium 
states, for which one has no a priori control of how close 
the system is to local equilibrium, has however remained 
difficult. One knows from approximate theories like fluc- 
tuating hydrodynamics [Q that such states exhibit long 
range correlations totally absent from equilibrium sys- 
tems (even at the critical point). These correlations have 
been measured experimentally in a fluid with a steady 
heat current B. Their derivation from a well defined 
macroscopic functional valid beyond local equilibrium is 
clearly an essential step in understanding pattern for- 
mation in more general stationary nonequilibrium states 
such as the Benard system. We report here what we be- 
lieve is the first exact derivation of such a functional for a 
nonequilibrium model which is relatively simple but ex- 
hibits the realistic feature of macroscopically long range 
correlations. 

For an equilibrium system, such as a lattice gas, in 
a unit cube containing L d sites with spacing 1/L (or a 
similar continuum system), at temperature T and chem- 
ical potential v, the probability of observing a speci- 
fied macroscopic density profile, with density p{x) at 
macroscopic position x in the unit cube, is given by 
P({p(x)}) ~ cxp[ - L d T cq ({p})]. Here T cq ({p}) = 
J \f(p(x)) — f(p(x))] dx, where the integration is over 
the unit cube, with f(p) the grand canonical free energy 
density and p the equilibrium density profile, obtained 
by minimizing J f(p(x))dx 0. The profile p(x) will be 
independent of x unless there is an external potential. 
We have suppressed the dependence of / and T eq on the 



constant temperature T, and assume that neither p nor 
p passes through a phase transition region at this tem- 
perature. 

In this letter we generalize the expression for T cq to the 
case of a system maintained, by contact with two bound- 
ary reservoirs at unequal chemical potentials i>o and v\, 
in a stationary nonequilibrium state with a constant par- 
ticle flux. We consider perhaps the simplest such system, 
the one-dimensional symmetric simple exclusion process 
on a lattice of N sites with open boundaries (To) . Each 
site i, i = 1, . . . , N, is either empty (r, = 0) or occu- 
pied by a single particle (n = 1). Each particle inde- 
pendently attempts to jump to its right neighboring site, 
and to its left neighboring site, in each case at rate 1. 
It succeeds if the target site is empty; otherwise noth- 
ing happens. At the boundary sites, 1 and N, parti- 
cles are added or removed: a particle is added to site 1, 
when the site is empty, at rate a, and removed, when 
the site is occupied, at rate 7; similarly particles are 
added to site N at rate 5 and removed at rate f3. We 
can think of sites i = and i = N + 1 as occupied with 
probabilities po — expfo/(l + exp^o) = a /(a + 7) and 
pi = exp 2/i/(l -I- expz^i) = 8/(13 + 6), independent of 
r. We will assume for definiteness that po > pi, but 
of course, due to the left-right symmetry, this is not a 
restriction. 

The probabilities of the microscopic configurations in 
the steady state may be calculated through the so-called 
matrix method jll]]. Here we ask for the probability of 
seeing a specified macroscopic density profile p(x), where 
< x < 1 and < p(x) < 1. By definition, this is the 
sum of the probabilities of all microscopic configurations 
consistent with the given profile. We shall not give a 
precise definition (see |||9)) of "consistent" here; roughly 
speaking we include all configurations r such that for any 
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y, z with < y < z < 

<5/v, with Sn — > as N — > 00. 

Our main result is a parametric formula for this prob- 
ability: for large N, P({p(x)}) - exp [-NT({p})}, with 
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H{p})= C dx{p{x)\o r l , ] 
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F{x) 



+ { l- p{x ))log[^^-\ + \og(I^\}. (1) 



1 - F{x) 



Pi - Po 



Here F is an auxiliary function determined by p(x): it is 
the decreasing solution of the differential equation 



p(x) = F(x) + 



F(x)[l-F(x)}F"{x) 



F'(x) 2 

satisfying the boundary conditions 

F(0) = po , F(l) = Pl 



(2) 



(3) 



It can be shown that such a solution exists and is (at 
least when p\ > and po < 1) unique |l5[ ]. 

We note that if one looks for a monotone function F, 
satisfying the constraint (|^), for which the expression 
Q({p},{F}) given by the right hand side of ([!]) is sta- 
tionary, one obtains (0) as an Euler-Lagrange equation: 



mip},{F}) 

SF(x) 



0. 



(4) 



We can in fact prove that this stationary point is a max- 
imum when < p\ < pa < 1 and expect this to be 
true in general: 



H{p}) = Bu P g({ P },{F}). 

F 



(5) 



Before explaining the main steps of the derivation of 
let us comment some of their consequences: 

(a) By its very nature T satisfies T{{p\) > 0, with 
equality only for p(x) = p{x), where p(x) = P q{1 — x) + 
p\x is the profile obtained with probability one in the 
limit N — ► co. Any other profile will have JF({p}) > 
and thus, for large N, exponentially small probability. 
When po = 1 or p\ = there are some profiles for which 
T = +oo; their probability is super-exponentially small 
in N. For examples, see (b) below and Jl5| . 

(b) For a constant profile p(x) — r, F satisfies F' = 
AF r (l - F) 1 ~ r , where A is fixed by (§), and 



H{P}) = log 
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dz 



Pi - Po 



(6) 



We see that T({p\) = oo if r = 1 and Pl = 0, or r = 
and po = 1. 

(c) Using M), one finds immediately that 



ST 
5p(x) 



log 



p(x) 1 - F(x) 
F(x) ' 1 - p(x) 



(7) 



This expression can then be used to find optimal profiles 
subject to various constraints. For example, setting (|t]) 
to zero implies that p(x) — F(x), leading to the most 
likely profile p(x) given in (a). 



(d) If we minimize T subject to the constraint of a 
fixed mean density J Q p(x)dx, the right hand side of (7) 
becomes an arbitrary constant, and together with (2) 
one obtains that the most likely profile is exponential: 
p(x) = A\ exp(#a;) + A2, the constants being determined 
by the value of the mean density and the boundary condi- 
tions ([}]). (This exponential form, which is the stationary 
solution of a diffusion equation with drift, was first sug- 
gested to us, for some special cases, by Errico Presutti). 

Similarly, if we impose a fixed mean density in k 
nonoverlapping intervals, with no other constraints, the 
optimal profile is exponential inside these intervals, linear 
outside, and in general not continuous at the end points 
of the intervals. 

(e) When the chemical potentials of the two reservoirs 
are equal, i.e., po = p\, the system is in true equilibrium 
with p(x) — p\. Eqs. (|l]-^|) have a well-defined limit for 
Pi / Po, with F(x) = po + (pi - po)x + (9((pi - p ) 2 ). 
It is also natural to consider a local equilibrium Gibbs 
measure corresponding to a spatially varying chemical 
potential |^|| which is adjusted to maintain the same 
optimal profile p(x). For this system the large deviation 
functional (free energy) is just J-" oq ({/o}), which has the 
explicit form 



^eq(M) 



p{x) \o\ 



P( x ) 

p{x) 



+ [l-p(»)]log (8) 

From (§: T({p}) > Q{{p},{p}) - ^ q (M), so that 
expressions (pl) and (0) are in general different 



H{p}) > F cq ({ P }), 



(9) 



with equality only for p{x) = p(x) or p = Pi] 

(f) Using the fact that the exponential is the optimal 
profile for the case of a fixed mean density in the entire 
interval, we may compute the distribution of M, the total 
number of particles in the system, in the steady state for 
large N. We find that the fluctuations of M predicted 
by (Q) are reduced in comparison to those in a system in 
local equilibrium (0) with the same p: 



hm iV- 1 [(M 2 ) SNS -(M) 2 ] 

- hm N-^iM'U-iMf]-^ 1 '^ 2 
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(10) 



We may also obtain (]I0|) by expanding p{x) about p(x) 
in (jf|). The result agrees with that obtained in jl4| di- 
rectly from the microscopic model and from fluctuating 
hydrodynamics 0. 

Derivation: let us now sketch the derivation of 
The probability of a configuration r = {n, . . . ,tjv} in 
the steady state of our model is given by 0| 



Pn(t) 



u>n{t) 



{W\(D + E) N \V) 



(11) 
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where the weights wjv(t) are given by 

LUN ( T ) = (w\nf =1 (nD + (l - n)E)\v) (12) 

and the matrices D and E and the vectors \V) and (W\ 
satisfy 

BE- ED = D + E , (13) 
(/3D-8E)\V) = \V) , (W\(aE-jD) = (W\ . (14) 

Although proving that &PJ) do give the weights in the 
steady state is rather easy fnj , there is not so far a simple 
physical interpretation of the matrices D and E or of the 
vectors \V) or (W\. 

To obtain the probability P Nlt ,,^ Nn (Mi, M 2 , M n ) 
that Mi particles are located on the first N% sites, M 2 
particles on the next N 2 sites, etc., we first calculate the 
sum f2jVi,...,Ar„(Mi, M 2 , . . . ,M n ) of the weights of all the 
corresponding configurations. The key to obtaining (|l|) 
is that the following generating function, which plays the 
role of the grand-canonical-pressure partition function, 
can be computed exactly: 



Z{\!,. 



, X n ; /zi, 
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&Ni,...,N n (Mi, M n ) 



(W\V) 



ol i n \ n D+ii n E 



(W\V) 



(15) 



where the sum is over all N, Mi > and the parameters 
Hi and Xi are conjugate to the Ni and Mi. To do the cal- 
culation we used repeatedly the following identity, which 
follows from (U3[), 



^xD+yE 



\xe y — ye x 
and the result is that 

a+b 



z = 



Po ~ Pi 
9 



exp 



(x - y)e x 
xe y — ye 1 



D 



(16) 



(17) 



where a ~ (a + 7) 1 , b = ((3 + 5) 1 , and 
g = -pi + Po e^i=i^ 



E 



1 



1 



(e 



A 3 ) 



(18) 



All the rest of the derivation consists in extracting from 
the exact expression (|l7]) the behavior of £1 and P fol- 
iar ge Ni's. 

First (0) gives the normalization factor in (p"l|), 



(w|y) 



i> + 6 + AQ 
r(a + 6)(p -pi) w 



(19) 



so that PjVi,...,iV„(-Mi, M2, . . . , M n ) is given simply by 
n Nu ... <Nn {Mi',M 2 ,...,M n ) divided by ©. 

Now the large Ni behavior of f2jy"i ....,w„ controls the 
position and the nature of the singularities of Z closest 
to the origin p,i = Xi = 0; conversely, this relation can 
be inverted |L5| to determine the asymptotic behavior 
of Ojvj ,...,Ar ra from the equation g = of the surface on 
which Z is singular, just as the growth of the coefficients 
of a power series of one variable is determined by the 
singularity nearest the origin. In particular, if one sets 
Nj = Ny 3 and M 3 = N 3 r 3 then, for large N and fixed 

Uit Tit 

logP j v 1 ,...,jv„(M 1 ,...,M Tt ) 
N 

a 

~ \og{p Q - p 1 )-^y j {\ogp, j /y j + r j \ogX j ) , (20) 



where y 3 and Tj are related to the parameters fix, . 
and Ai, . . . , A„ by 



dg 



EI 



dg 



'log Mi 



<9 log A j 

a log fi j 



1 Mr. 



(21) 



with all derivatives calculated on the manifold g = 0. 

Equations ( pp| ) and ( pl[ ) determine J 7 in a parametric 
form. As the (i/s and the Aj's vary, they give the sizes 
of the boxes yj , their particle densities r 3 , and the corre- 
sponding probabilities. Note that the parameters a and 
b do not appear in ( |l8| ) and therefore in the equation 
g = 0; so for large N, only po and p\ remain relevant. 

This parametric form can be simplified by replacing 
the role of the fii and A^ by a single sequence of pa- 
rameters Gi. Let us define the constant C by C = 
E™=i 9g/dlogpi and the sequence Gi by 



Gi = C 



G 



n+l 



c~ 



(22) 



It follows from ( p2| ) that = log(Gi/Gi+i)/(l — A^) and 

1 1 y l 



X t - 1 G l+1 log(G t /G l+1 ) 



Po 



±(--—y 



Gj G J+1 '\og{G 3 /G J+1 ) 
The condition that g = becomes 



(23) 



^ 1 I 



Vi 



^G 3 G 3+1 >log(G 3 /G J+1 )> 



(24) 



which can be thought of as an equation which determines 
G n +i in terms of G\, . . . , G n . 

Once the Ai's are known, one gets for the n's and the 
large deviation function 
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A, 



A, 



d — G. 



i+l 



1 - Aj 

and 

log[P Nl ,..., Nn (M u . 



(1-A,) 2 
>M n )] 



(25) 



N 



(26) 



n 

jri i- Lj/j(po - pi) 

-y 2 log(l - Aj) + y,r i log(A. i )}- 

To take the limit N — > oo followed by n — ► oo, yi — > 0, 
we introduce a continuous variable < x < 1, and 
let Xi = iji + y 2 + . . . + jji- All the discrete sequences 
can now be thought of as functions of xi,...,x n with 
Gi = G(xi), Xi = X(xi) and r 4 = p(xi), i = 1, ...,n. 
Then, extending G, A, and p to be smooth functions of 
x, so that d+i — d ~ yiG'(xi), one finds that (HHJ) 
become 



1 



-1 



A(x) - 1 
Pi 



G'(z) 
Po - 



Po 



du 

W) 



du 



(27) 
(28) 



At this stage it is more convenient to replace G{x) by 
the function F(x) defined by F(x) = p + L G(u) _1 du. 
The expression (|27]) for A (a;) becomes then 



X(x) 



F'(x) 2 
F"{x) 



F{x). 



(29) 



Using the above relations we may rewrite (J25|) and (|26 
in terms of F 7 obtaining (|l]) and (||); the boundary condi 
tions (^) come from ((2^). The monotonicity of F follows 
from the uniqueness of the sign of G (see (|22])). 

Possible extensions: It would be nice to give a physical 
interpretation of the Euler-Lagrange equations di- 
rectly from a macroscopic or semimacroscopic approach. 
Such an approach might allow an extension to the large 
deviation function for time-dependent profiles and per- 
mit one to unify various results recently obtained for the 
large deviation functions of time dependent profiles or 
currents §@[[|- 

The matrix method (with modified rules (|l3T4)) also 
applies pd| ] to the asymmetric simple exclusion process, 
in which particles jump to the left and right at different 
rates. We are at present trying to generalize our cal- 
culations to this case. We have already calculated the 
probability of a specified mean density p (total number 
of particles) for the special case in which particles jump 
only to their right and with only input, with rate a > 1, 
at the left and output, with rate /3 > 1, at the right: 

N- 1 log P(p) ~ -2[ /9 log2p + (1 - p) log 2(1 - p)\ . (30) 

For a Bernoulli measure at uniform density 1/2, the result 
would be the same except for the factor 2 in front of 



the whole expression, which here again expresses the fact 
that fluctuations are reduced by long range correlations. 

Conclusion: For the simple model we studied here, 
the large deviation function T (which extends the notion 
of free energy to nonequilibrium systems) is a nonlocal 
functional (|l|-§) of the density p. This implies that the 
probability of a given profile is not the product of proba- 
bilities for, say, the left and right halves of the system, in 
contrast to the situation for equilibrium systems. Since 
long range correlations are expected for general station- 
ary nonequilibrium states and have even been measured 
experimentally ||, the nonlocal nature of F({p}) is pre- 
sumably also general, in contrast with the conjecture 
that T is a local function of the hydrodynamic variables. 

We thank E. Presutti for very helpful discussions. 
The work of J. L. Lebowitz was supported by NSF 
Grant DMR-9813268, AFOSR Grant F49620/0154, DI- 
MACS and its supporting agencies, and NATO Grant 
PST.CLG.976552. J.L.L. and E. R. Speer acknowledge 
the hospitality of the I.H.E.S. in the spring of 2000, where 
this work was begun. 



[10 
[11 



[12] 

[13] 
[14] 
[15] 



L. Onsager, S. Machlup, Phys.Rev. 91 1505 (1053) ; 91 
1512 (1953) 

C. Kipnis, S. Olla, and S. R. S. Varadhan, Commun. Pure 
Appl. Math. 42, 115 (1989). 

C. Kipnis and C. Landim, Scaling Limits of Interacting 
Particle Systems (Springer- Verlag, Berlin, 1999). 
R. Graham, in Order and Fluctuations in Equilibrium 
and Nonequilibrium Statistical Mechanics, ed. G. Nicolis, 
G. Dewel, and J. W. Turner (Wiley, New York, 1981). 

G. Eyink, J. Stat. Phys. 61, 533 (1990). 

L. Bertini et ai, Phys. Rev. Lett. 87, 040601 (2001) 
J. R. Dorfman, T. R. Kirkpatrick, and J. V. Sengers, 
Annu. Rev. Phys. Chem. 45, 213 (1994). 
W. B. Li et al, Phys. Rev. Lett. 81, 5580 (1998). 

H. Spohn, Large Scale Dynamics of Interacting Particles 
(Springer- Verlag, Berlin, 1991). 

G. Eyink, J. L. Lebowitz, and H. Spohn, CMP 132, 253 
(1990); 140, 119 (1991). 

B. Derrida, M. R. Evans, V. Hakim and V. Pasquier, 
J. Phys. A 26, 1493 (1993); F. H. L. Essler, V. Rittenberg 
J. Phys. A 29, 3375 (1996); N. Rajewsky, M. Schrecken- 
berg Physica A 245, 139 (1997); T. Sasamoto, J. Phys. A 
32, 7109 (1999); R. A. Blythe, M. R. Evans, F. Colaiori, 
F. H. L. Essler, J. Phys. A 33, 2313 (2000). 
B. Derrida and J. L. Lebowitz, Phys. Rev. Lett. 80, 209 
(1998). 

M. Prahofer, H. Spohn, Phys. Rev. Lett. 84, 4882 (2000). 

H. Spohn, J. Phys A. 16, 4275 (1983). 

B. Derrida, J. L. Lebowitz, and E. R. Speer, in prepara- 
tion. 



4 



